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Abstract 

In this paper, we propose a numerical method to solve isotropic elliptic equations on point 
cloud by generalizing the point integral method. The idea of the point integral method is to 
approximate the differential operators by integral operators and discretize the corresponding 
integral equation on point cloud. The key step is to get the integral approximation. In this 
paper, with proper kernel function, we get an integral approximation for the elliptic operators 
with isotropic coefficients. Moreover, the integral approximation has been proved to keep the 
coercivity of the original elliptic operator. The convergence of the point integral method is also 
proved. 


1 Introduction 

Nowadays, data plays more and more important roles in science and engineering. In many problems, 
data is usually represented as a collection of points embedding in a high dimensional Euclidean 
space. Processing and analysis of the point cloud data is essential in many applications, such as 
machine learning mm and image processing pH .'50 . 

In many applications, the point cloud data lies in a manifold whose dimension is much lower 
than the ambient Euclidean space. The low dimensionality is an important feature we could exploit 
to analyze the data. One example is the low dimensional manifold model (LDMM) in image 
processing m- In this model, the original image is cutted to many overlap patches. The collection 
of all patches consists of a point cloud in Euclidean space. It is found that for many natural 
images, the patch set usually samples a low dimensional manifold which is called patch manifold. 
The dimension of the patch manifold is used as a regularization to processing the image. Based on 
differential geometry and variational method, this model is reduced to solve Laplace equation on 
patch set. The key point in LDMM is to solve this Laplace equation accurately and efficiently. 

Beside the data analysis, solving PDEs on manifold also appears in many physical problems, 
such as material science PE], fluid flow 119. i2T] . biology and biophysics [21 US, 29, [2]. To solve 
PDEs on manifold, many methods have been developed, especially on 2D surfaces, including surface 
finite element method [16], level set method mm, grid based particle method [25, 23] and closest 
point method [321 [28 ]. However, these methods need extra information besides the point cloud, for 
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instance, meshes, level set function and closest point function. These information is not easy to 
obtain from point cloud when the dimension of the manifold is high. 

Recently, Liang et al. proposed to discretize the differential operators on point cloud by local 
least square approximations of the manifold EZ3- Their method can achieve high order accuracy 
and enjoy more flexibility since no mesh is needed. In principle, it can be applied to manifolds with 
arbitrary dimensions and co-dimensions with or without boundary. However, if the dimension of 
the manifold is high, this method may not be stable since high order polynomial is used to fit the 
data. Later, Lai et al. proposed local mesh method to approximate the differential operators on 
point cloud [23]. The main idea is to construct mesh locally around each point by using K nearest 
neighbors. The local mesh is easier to construct than global mesh. Based on the local mesh, it 
is easy to discretize differential operators and compute integrals. However, when the dimension of 
the manifold is high, even local mesh is not easy to construct. 

The original point integral method for Laplace equation is closely related with the graph Lapla- 
cian mm- Graph Laplacian has been widely used in many problems. It is observed in [5] 122] 20, 35] 
that the graph Laplacian with the Gaussian weights well approximates the Laplace-Beltrami oper¬ 
ator when the vertices of the graph are assumed to sample the underlying manifold. When there 
is no boundary, Belkin and Niyogi [6j showed the spectra of the graph Laplacian with Gaussian 
weights converges to that of Laplace-Beltrami operator. Recently, Singer and Wu [38] showed the 
spectral convergence of the graph Laplacian in the presence of the Neumann boundary. 

Inspired by the graph Laplacian and the nonlocal diffusion, we developed the point integral 
method for Poisson equation on point cloud [26], [33, 34] • 

-A M u(x) = f(x), x E M, 


where = div(V) is the Laplace-Beltrami operator in AL 

We assume that Af E C°° is a compact ^-dimensional manifold isometrically embedded in 
with the standard Euclidean metric and k < d. If Ai has boundary, the boundary, <9Af is also a 
C°° smooth manifold. 

Let <h : Q C —> M C be a local parametrization of M and 9 E fh For any differentiable 

function / : J\4 —> M, define the gradient on the manifold 


vnm) = E 9 ij (o)^-(6) df ^} 0)) (9), 


i,j =1 


dOi 


d6i 


( 1 . 1 ) 


and for vector field F : M. —> T x M. on AL where T X M. is the tangent space of M. at x E AL the 
divergence is defined as 


div(F) 


d 


^atu% ds - 


EE» VdoG^FHm) 


Wi 


( 1 . 2 ) 


where ,k = G , detG is the determinant of matrix G and G{9) = (gij)i,j=i,- ,k is the 

first fundamental form which is defined by 


d 


9ij(0) ~ El 

k =1 


d<s> k d<s> k 


i,3 = 1,' 


, m. 


(1.3) 


and (E 1 (x), • • • ,F d (x)y is the representation of F in the embedding coordinates. 
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The main idea of the point integral method is to approximate the Poisson equation by the 
following integral equation: 


- / A 7 v[ti(y)i2 i (x,y)d/xj 






f — du 

i? i ( x ,y)(n(x) - n(y))d/i y - 2 / R t (x,y) — (y)dr y , 


IdM 


dn 


where n is the out normal of A4. i?i(x, y) and Rt(x,y) are kernel functions given as follows 


.R t (x,y) = C t R 


|x- y| 
At 


, R t (x 1 y) = C t R 


|x - y| 
4 1 


(1.4) 


where C t = 
over [0,+oo). And 


is 


the normalizing factor. R E C' 2 (M + ) be a positive function which is integrable 


R(r) = 


P+OO 


R(s)ds. 


There is not any derivatives in the integral equation. It is easy to be discretized from point clouds 
using some quadrature rule. In [33, [33], we proved the convergence of the point integral method 
for Poisson equation with Neumann and Dirichlet boundary condition. 

In the point integral method, we only need the point cloud to discretize the differential operator. 
This gives PIM great flexibility to fit the requirements in variety of applications. However, one 
limitation of the point integral method is that it only applies on Laplace-Beltrami operator. In many 
problems, we need to discretize other differential operators besides Laplace-Beltrami operator. In 
this paper, we generalize the point integral method to isotropic elliptic operators. Isotropic elliptic 
operators are also widely used in many problems. One example is the nonlocal total variation 
minimization on point cloud, in which we need to solve an optimization problem, 


min ||Vu||ii(_A 4 ), subject to: T(u) = b 


where V is the gradient in A4, T is the measurement operator related with the application, b is the 
observation and 

l|Vu|| L i( M) = [ |Vu(x)|dx 

J M 

Using standard variational approach, the solution of above optimization problem can be given by 
solving a nonlinear elliptic equation, 


—div 


( V ^(x) \ 

V|Vu(x)|y 


/(x). 


where /(x) is a known function. Apparently, this equation can be solved by solving a sequence of 
isotropic elliptic equation iteratively. 

In this paper, we consider to solve elliptic equations with isotropic coefficients on manifold M, 


—div(p 2 (x)Vu(x)) = /(x), x E M 


(1.5) 


The coeffcients p(x) and source term /(x) are known smooth functions of spatial variables, 


i.e. 


peC\M), feC\M). 


The elliptic condition makes that there exist generic constants cq, ci >0 such that for any x E A4, 


0 < cq < p(x) < Cl < oo, 
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The key observation in this paper is the integral approximation of isotropic elliptic operators given 
as following 


[ div(p 2 (y)Vu(y)) ^ X ^ rf/i y [ 

Jm p(y) t Jm 


R t {-x, y)(u(x) - u{y))p{y)dp y (1.6) 



where the kernel functions Rt and Rt are same as those in (1.4). The main advantage of this 


integral approximation is that there is no differential operator inside. Using this approximation, we 
transfer the numerical differential to numerical integral which is much easier to compute on point 
cloud. Based on this integral approximation, we are able to develop the point integral method to 
isotropic elliptic equations. 

Similar integral approximation is also widely used in nonlocal diffusion and peridynamic model 
cam isnm Eg. The integral approximation is easy to implement on point cloud, since it has no 
derivatives inside. Moreover, the point integral method also has very good theoretical property. It 
is proved that the coercivity of the original elliptic operator is partially preserved and this partial 
coercivity implies the convergence of the point integral method. 

The rest of the paper is organized as following. In Section 2, we introduce the point integral 
method for isotropic elliptic operator with Neumann and Dirichlet boundary condition. The con¬ 
vergence analysis is given in Section 3. Several numerical examples are presented in Section 4. The 
conclusion remarks are made in Section 5. 


2 Point Integral Method for Isotropic Elliptic Equations 


In this section, we introduce a numerical method for isotropic elliptic equation on point cloud based 


on the integral approximation (1.6). 


To simplify the notation, we introduce an integral operator, 



( 2 . 1 ) 


where Rt is the kernel function given in (1.4). 


2.1 Neumann Boundary 

First, we consider the Neumann problem, 



( 2 . 2 ) 


Using the integral approximation (1.6), the solution of the Neumann problem (2.2) can be obtained 


approximately by solving an integral equation 



(2.3) 


with id. 

The eigenvalue problem is also solved by a generalized eigenvalue problem 



(2.4) 
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2.2 Dirichlet Boundary 


The Dirichlet problem is more involved in point integral method, since the normal derivative, 
is not known. 

f —div(p 2 (x)Vu(x)) = /(x), xGM, . 

\ u(x) = g(x), x e dM. 

Here, we use the same idea as that in [26] to deal with the Dirichlet boundary. 


2.2.1 Robin Approximation 


The simplest way is using Robin boundary to approximate the Dirichlet boundary. More specifically, 
we consider the following Robin problem 


—div(p 2 (x)Vu(x)) = /(x), x e M, 


u(x) + /3^(x) = g(x), xedM. 


( 2 . 6 ) 


where 0 < f3 <C 1 is a small parameter. It is easy to show that as (3 —> 0, the solution of the Robin 
problem, (2.6), converges to the solution of the Dirichlet problem, (2.5). 

For the Robin problem, the integral approximation (1.6) is applicable to give an integral equa¬ 
tion, 


L ^( x )-| [ (g(y) -u(y))Rt(x,y)p(y)dT y = [ ^-R t (x, y)dg y . (2.7) 

P JdM Jm p{ y) 

Similarly, we also get an approximation of the eigenvalue problem, 

L t:u(x) + w / u(y)R t (x,y)p(y)dr y = A f u(y) ^^^ dp y (2.8) 

P Jom Jm P(y) 


2.2.2 Iterative Solver based on Augmented Lagrangian Multiplier 


In the Robin approximation, the parameter /3 has to be small to get good approximation, while the 
linear system becomes ill-conditioned. To alleviate this difficulty, we could use an iterative method 
based on the Augmented Lagrange method (ALM). 

It is well known that the Dirichlet problem can be reformulated to be following constrained 
optimization problem: 

min i [ p 2 (y)\Vv{y)\ 2 dn y + [ /(y)u(y)d/r y , (2.9) 

subject to: u(x)| aA ^ = g(x), 


Applying the ALM method to the problem (2.9), we get an iterative method, in each step, an 
unconstrained optimization problem is solved, 


min -f p 2 ( y)|Vu(y)| 2 d/i y + / f(y) ■ v(y)dg y 

+ [ w k (y) • (g(y) - v(y))p 2 (y)dr y + f (ff(y) - u(y)) V(y)dr y . (2.10) 
JdM Z P J dM 

Using the variational method, one can show that the solution to ( |2.10[ ) is exactly the solution to 
the following Robin problem: 


f div(p 2 (x)Vv(x)) = /(x), xgM, . 

\ v(x) + P§^(x) = g(x) + f3w k (x), x e dM. J 

This Robin problem is solved by the integral equation. Notice that, the parameter (3 is not neces¬ 
sarily small. Usually, we set (3 = 1 . Thus, the linear system is not ill-conditioned. 
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Procedure 1 ALM for Dirichlet Problem _ 

1 : k = 0 , w° = 0 . 

2: repeat 

3: Solving the following integral equation to get v k , 

L t v k { y) ~ i [ (g(y) - v k {y) + /3w k (y))R t (x,y)p(y)dTy = f ^^(x, y)d/x y . 

p JdM Jm p{ y) 

4 : W k+1 = w k + jj(g — ( v k \dM )), k = k + 1 

5: until ||g - (n fc_1 |aM)|| = 0 

6 : U = V k 


2.3 Discretization 


The main advantage of the integral equations is that they are easy to discretize over the point cloud 
since there is not any derivatives inside. 

Assume we are given a set of sample points P sampling the submanifold Ai and a subset S C P 
sampling the boundary of Ai. List the points in P respectively S in a fixed order P = (xi, • ■ ■ , x n ) 
where x* E M d , 1 < i < n, respectively S = (si,--- ,s nb ) where sj E P. In addition, assume 
we are also given two vectors V = (Vi, 

A = (Ai, ■ 


,V n Y where Vi is an volume weight of x,; in Ai, and 
A nb Y where Ai is an area weight of s* in dAi. In this point cloud data (P, S, V, A), 


the integral equation (2.3) can be disretized as 


1 


X! k\(x,.xj)(ui - Uj)pjVj = 2 Rt(xi,Bj)bjPjAj + !>’/(x,. xj)fjVj/pj 


( 2 . 12 ) 


x,eP 


Si GlS" 


x,eP 


where pj = p(xj), fj = f(xj), b t = b(s t ), j = 1, ■ ■ • , |P|, l = 1, ■ ■ ■ , |5|. 

The other integral equations and corresponding eigenvalue problems can be discretized conse¬ 
quently. 


Remark 2.1. The integral approximation (1.6) also holds if the parameter t depends on x, i.e. 

|2 ' 


‘ L J R 


|x - y| 

4t(x) 


(«(x) - u(y))p(y)dp y (2.13) 


-2 f p(y)R 

JdM 


|X - y| 

4 t(x) 


P(y)dr y . 


Based on above approximation, in the computation, we can choose t adaptive to the distribution of 
the points. 


3 Convergence Analysis 

In this section, we analyze the convergence of the point integral method for isotropic elliptic equa¬ 
tion. To make the theoretical analysis concise, we only consider the homogeneous Neumann bound¬ 
ary conditions, 


f -div(p 2 (x)Vri(x)) = /(x), x E A4, 
{ g*(x) = 0, xedAi. 


(3.1) 
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The corresponding numerical scheme is 



(3.2) 


The analysis can be easily generalized to the non-homogeneous boundary conditions. The conver¬ 
gence of Dirichlet problem can be proved also following the similar procedure as that in [31]. 

3.1 Main Result 

We will prove that the solution given by the point integral method converges to the exact solution 
as the point cloud (P, V) converges to the manifold At. Before giving the result of the convergence, 
we need to clarify the meaning of the convergence of the point cloud (P, V) to the manifold At. 

First, we introduce an index to measure the distance between the point cloud (P, V) and the 
manifold At, which is called integral accuracy index , denoted as h(P, V,At). 

Definition 3.1 (Integral Accuracy Index). For the point cloud (P, V) which samples the manifold 
At, the integral accuracy index h(P, V,At) is defined as 




feCHM) |supp(/)|||/||ci(A^) 


where ||/||c 1 (A'i) = ll/lloo + ||V/||oo and |supp(/)| is the volume of the support of f. 

Using the definition of integrable index, we say that the point cloud (P, V) converges to the 
manifold At if h(P, V, At) —> 0. In the convergence analysis, we consider the case that /i(P, V, At) 
is small enough. 

Remark 3.1. In some sense, h(P, V, At) is a measure of the density of the point cloud. If the 
point cloud is uniformly distributed on the manifold, from central limit theorem, h(P, V,At) ~ 
0(1/ \/\P\) where |P| is the number of point in P. 

Remark 3.2. To consider the non-homogeneous Neumann boundary condition or Dirichlet bound¬ 
ary condition, we have to also assume that h(S, A,dM.) —> 0, where S is the point set sample the 
boundary <9At and A is the corresponding volume weight on the boundary <9At . 

To get the convergence, we also need some assumptions on the regularity of the submanifold 
At and the integral kernel function R. 

Assumption 3.1. • Smoothness of the manifold: At,<9Al are both compact and C°° smooth 

fc-dimensional submanifolds isometrically embedded in a Euclidean space M d . 


• Ellipticity: there exist generic constants co,ci > 0, such that Co < p(x) < c\ and p(x) E 

C' 1 (At). 


• Assumptions on the kernel function P(r): 

(a) Smoothness: R E C 2 (M + ); 

(b) Nonnegativity: R(r) > 0 for any r > 0. 

(c) Compact support: R(r) = 0 for Vr > 1; 

(d) Nondegeneracy: > 0 so that P(r) > <5q for 0 < r < ),. 
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Remark 3.3. The assumption on the kernel function is very mild. The compact support assumption 
can be relaxed to exponentially decay, like Gaussian kernel. In the nondegeneracy assumption, 1/2 
may be replaced by a positive number do with 0 < do < 1. Similar assumptions on the kernel 
function is also used in analysis the nonlocal diffusion problem Hi- 


All the convergence analysis in this paper is based on above assumptions. In the statement of 
the theorems, above assumptions are omitted to make the statements more concise. 

The other issue we have to address is that how to compute the difference between the discrete 


solution and the analytic solution. The solution of the discrete system (|3.2|) is a vector u defined on 
P while the solution of tl 
for any solution u = (ui, ■ 


P while the solution of the problem (3.1) is a function defined on Ai. To make them comparable, 

\t 


) 'U' T 


n = \P\ to the problem (3.2), we construct a function on A i 


/ f (u)(x) = 


_ Ex jg P R t ( x > Xj)"jP.i V j + * Exjgp X ilf., V i.P.i 


V.. 




■ S.nV 


(3.3) 


It is easy to verify that If(u) interpolates u at the sample points P, i.e., If(u)(xj) = uj for any 
x.j E P. The following theorem guarantees the convergence of the point integral method. 


Theorem 3.1. Let u be the solution to problem (3.1) with f E C' 1 (Ad) and let the vector u be the 
solution to the problem (3.2). Then there exists constants C and To depend on M. andp(x), such 
that for any t, < To, 


u — I f (u)\\ H i {M) < C ^t 1//2 H—ll/llcR-M)- ( 3 -4) 


where h(P, V,A1) is the integral accuracy index. 


3.2 Proof of Convergence 

Roughly, the proof the convergence includes two parts: estimate of the truncation error Lt{u— If (u)) 
and the stability of the integral operator Lt. Here Lt is the integral operator in (2.1) , u(x) is the 
solution of the problem (3.1) and u is the solution of the problem (3.2). 

This strategy is standard in numerical analysis. It is well known that consistency together 
with stability imply convergence. On the other hand, the point integral method has some special 
structures both in truncation error and stability, which makes the analysis a little more involved. 
First, we have following theorem regarding the stability of the operator Lt. 

Theorem 3.2. Let u(x) solves the integral equation 


L t u = r(x) 

where r E H 1 (A4) with j M r(x)p(x)dp x = 0. Then, there exist constants C > 0, To > 0 independent 
on t, such that 

< c (IMI L 2 (M) + i ll^ 7r ll L 2 (M)) 

as long as t <Tq. 

To use above stability result, we need L 2 estimate of L t (u — If (u)) and \7L t (u — If(u)). In the 
analysis, we split the truncation error L t {u — If(u)) to two terms, 


L t (u - If(u)) = L t (u - u t )) + L t (u t - If (u)) 













where ut is the solution of the integral equation 


t 


'M 


i?t(x,y)(u(x) - u(y))p(y)d/i y = 


f t , ^Rt(x,y) 

/ / y - 

/A4 ?y 


(3.5) 


For the second term, we have following estimate. 


Theorem 3.3. Let rq(x) be the solution of the problem (3.5) and u be the solution of the prob¬ 
lem (3.2). If f E C l (M) , f/ien there exists constants C,Tq depending only on M. and the coefficient 
p(x), so that 


\\Lt ( 7 fU — ut) \\l 2 (M) < ^ £3/2 ’--ll/llch^i)’ 

|VL t (JfU — u,*) ||l 2 (X) < ^ ' t 2 - hfWc^M)- 

h(P,~V ,M) 


(3.6) 

(3.7) 


as long as t <Tq and 1 —- < To, h(P, V, Af) is the integral difference index in Definition 


3.1 


The error term Lt(u — ut )) is a little more complicated. It has two parts, one is the interior 
term and the other is the boundary term. We need to estimate these two terms separately to get 
better estimation of the convergence rate. 


Theorem 3.4. Let u(x) be the solution of the problem (3.1) and ut(x.) be the solution of the 
corresponding integral equation (3.5). Let 


d .. 

hd = V / n J (y)(x - y) • V(V J u(y))Ti(x,y)p(y)dr y , (3.8) 

JdM 


and 


Lt(u Ut) — I in T Ibd- 


where n(y) = (n 1 (y), • • • ,n d ( y)) is the out normal vector of dM. at y, V J is the jth component of 
gradient V. 

If u G H 3 (. M) , then there exists constants C,Tq depending only on M andp(x.), so that, 

\\Iin\\L 2 (M) — ^ 1 ^ 2 \\ u \\h 3 {M)i \\VIin\\L 2 (M) — C\\u\\h3(M)i ( 3 - 9 ) 

as long as t < To. 

Using the definition of the boundary term Im, ( |3.8[ ), it is easy to check that 
\\hd\\ L 2 {M) = 0(i 1/4 ), \Nhd\\ L 2 (M) = or 1 / 2 ), 

Based on this estimation, Theorem |3.2| and Theorem |3.4| give that 

\W~u t \\m(M) = 0{t 1/A ). 


This proves the convergence, however the convergence rate is relatively low. This low rate comes 
from the boundary term. From interior term only, the rate is \/i. Notice that the boundary term 
has a specific integral formula given in (3.8). Using this formula, we know that the boundary term 
concentrates in a small layer adjacent to the boundary whose width is of the order of \ft and vanish 
in the interior region. Utilizing this special structure, we could get better convergence rate with 


the help of a stability estimate specifically for the boundary term, which is given in Theorem 3.5 
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Theorem 3.5. Let u(x) solves the integral equation 


Ltu = / b(y) • (x - y).Rt(x, y)p(y)dr y - b 
JdM 


where \M\ P = f M p(x)dp, x and 


b = 


b(y) • (x - y)i?t(x, y)j?(y)dr y dx 


\M\ p \JdM 

Then, there exist constant C > 0, To > 0 independent on t, such that 

IMI H^iM) — l|b||^i(>i)- 

as long as t <Tq. 


Based on above four theorems, it is easy to prove Theorem 3.1 Using Theorem 3.3 and Theorem 
3.2[ we get 

\\ u t ~ ~ O I -— 


Applying Theorem 3.2 to the interior term in Theorem 3.4 and Theorem 3.5 to the boundary term 
respectively, we have 

\\ u ~ = O 


Putting above two inequality together, Theorem |3.1| is proved. 
Next, we prove Theorem 3.2 , [3731 [3 d] and |3.5| respectively. 

3.3 Proof of Theorem 13.41 


Let r(x) = L t u — L t ut where u and ut are the solution of (3.1) and (3.5) respectively. Using 
integration by parts, we have 


r x =- 


T [ ^t(x,y)(u(x) - u(y))p(y)dn y - f div(p 2 (y)\/u{y)) Rt ^'^ d/x y (3.10) 

1 jm jm p\ y) 


du 


-2 R t (x,y) — (y)p(y)dT y 
JdM 011 


7 [ («( x ) - u{y) - (x - y) • Vu(y))R t (x, y)p(y)dp y 
1 Jm 


A M u(y)R t {x, y)p(y)dn y 


Jm 

The main idea of the proof is the Taylor expansion, 

u(x) - u( y) - (x - y) • V«(y) = ^(x - y) T ■ H u (y) • (x - y) + 0(|x - y| 3 ) 
where H u (y) is the Hessian matrix of u at y. 

Using the integration by parts, the second order term actually gives Laplace-Beltrami operator 
which cancel with the second term in (3.10). 

In manifold, the Taylor expansion and integration by parts are more complicated. To make 
the whole idea rigorous, we need to introduce a special parametrization of the manifold J\4. This 
parametrization is based on following proposition. 
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Proposition 3.1. Assume both M and dM are C 2 smooth and a is the minimum of the reaches 
of M and dM. For any point x G M, there is a neighborhood U C M of x, so that there is a 
parametrization $ : 12 C — > U satisfying the following conditions. For any p < 0.1, 

(i) 12 is convex and contains at least half of the ball -B$-i( x )({?cr), i.e., vol(Q n (s' 7 )) > 

\{^(j) k Wk where Wk is the volume of unit ball in M fc ; 

(%%) B x (^a)UM C U. 

(iii) The determinant the Jacobian of & is bounded: (1 — 2 p) k < < (1 + 2 p) k over 12. 

(iv) For any points y, z G U, 1-2 p< |$ „ 1( j^_ 1(z)| <1 + 3 p. 


This proposition basically says there exists a local parametrization of small distortion if (M, dM) 
satisfies certain smoothness, and moreover, the parameter domain is convex and big enough. The 
proof of this proposition can be found in [33j and for the sake of completeness, we give the proof 
in the supplementary material. Next, we introduce a special parametrization of the manifold M. 

Let p = 0.1, a be the minimum of the reaches of M and dM and d = pa/ 20. For any x G M, 
denote 

Bi = {y G M : |x — y| < 5} , M x = {y G M : |x - y| 2 < At] (3.11) 

and we assume t is small enough such that 2 y/t < 5. 

Since the manifold M is compact, there exists a <5-net, Ms = {q* G M, i = 1, ■ ■ • ,N}, such 
that 

N 

Mc[jBl 

1=1 

and there exists a partition of .Ad, {Oi, i = 1, • • • ,N}, such that 0 % n Oj = 0, i ^ j and 

N 

M = \JO u Oi C Sq., i = l,...,N. 

2—1 

there exist a parametrization <f>j : 12* C —> Ui C M, i = 1, • • • , N, 

such that 


Using Proposition 


3.1 


1. (Convexity) C Ui and 12* is convex. 

2. (Smoothness) <1+ G C 3 (12j); 

3. (Locally small deformation) For any points 0\,0 2 £ 1 

\ \0i - e 2 1 < 11 ^( 0 ,) - ^(0 2 )|| < 2 ift - e 2 


Using the partition, {Oi, i = 1, ■ • • ,N}, for any y G M, there exists unique J(y) G {1, • • • ,N}, 
such that 

y G °J(y) C < (y) - (3.12) 

Moerover, using the condition, 2 y/t < 5, we have My C C Uj( y y Then <hj^(x) and 

(y) are both well defined for any x G M y . 
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(3.13) 


Now, we define an auxiliary function, jy(x, y) for any y E M, x E M y . Let 


i.e. 


?(x,y) = $ J( y) (x) - $ J(y) (y) € K , r/(x,y) = f(x,y) • 5<L J(y) (a(x, y)) E R , 

where a(x,y) = <&J^(y) and d is the gradient operator in the parameter space, 

= een,cR> 

Now we state the proof of Theorem |3.4| 


Proof. First, we split the residual r(x) in ( 3. 10[ ) to four terms 

r(x) =ri(x) + r 2 (x) + r 3 (x) - r 4 (x) 

where 


r i{x) = 


'M 


u{x) - u(y) - (x - y) • Vu( y) - (V l V J u(y)) ) i2*(x, y)p(y)dp y , 


^ 2 ( x ) = ^ / ?7*77' 7 (V*V J u(y))i?j(x, y)p(y)d/r y — / r ? *(V l V%(y)V J i?i(x, y)p(y)d^ y , 


/A4 


/X 


= 


/X 


? 7 *(V*V- 7 tt(y)V J i?i(x,y)p(y)d/r y + / div (r/*(V*Vu(y)) i? t (x,y)p(y)d/x ; 


/X 


y> 


r 4 (x) = [ div (77*(V*Vu(y)) i? i (x,y)p(y)d/x y + [ A M u(y)R t (x,y)p(y)dp y . 

JM JM 

where V*, i = 1, ■ • • , d is the ith component of the gradient V, rf, i = 1, • • • , d is the ith component 
of r](x, y) defined in (3.13). To simplify the notation, we drop the variable (x, y) in the function 

d( x ,y)- 

Next, we will prove the theorem by estimating above four terms one by one. First, we consider 
r\ . Let 

d(x, y) = u(x) - u(y) - (x - y) • Vu(y) - ^ i ^(V l V J u(y)). 


we have 


/ l r i( x )| 2 d/i x = / / R t (x,y)d(x,y)p(y)dp y 

JM JM JM 


>M JM 
< (rnaxp(y)) 

y jm \jm 


dp y 


< C 


im Jm 


Rt{x,y)dp, y S j d? 4 (x,y)|d(x,y)| 2 d/i y ^ dp* 
R t {x,y)\d{x,y)\ 2 dp y dp x 


and 


/ / Rt(x,y)\d(x, y)\ 2 dp y dp x = V / / ^t(x, y)|d(x, y)| 2 d/z y d/i x 

Jm Jm “7 ix JOi 

= I ( [ R t(^y)\d{x,y)\ 2 dp^\ dp y . 

i=l Jo ‘ \ J M* y ) 
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Using Newton-Leibniz formula, we get 


d(x,y) = «(x)-u(y)-(x-y)-Vu(y)--r ? y(V*VUi(y)) 

= CC' 1 j J J + S3Si£)<9j/$ J, (a + S 3 S 2 si^)V j 'V j u(^(a + s 3 s 2 si£))) ds 3 ds 2 dsi 

= cc'c" f [ [ sls 2 di&(a + s 3 si£)di"di'&'(a + s 3 s 2 si^)'V j '’V j u(^(a + s 3 S2Si^))ds 3 ds2dsi 
J 0 J 0 J 0 

+cc'c" [ [ [ sld i »di$ j (a + S3Sit,)di’&\a + S3S2Si£,)\7 j 'V j u($(a + s 3 S2Si£))ds 3 ds 2 dsi 
Jo Jo Jo 

»i />i r i 

eld'd" I 2 


+CCC I / / + S 3 S2SlO^'^ (« + S 3 S2SlO^"^ (a + S 3 S 2 Sl£) 

Jo JO Jo 

V J V J V%(<l>(a + S3S2SiO)ds3ds2dsi 

Here, d> 1 , z = 1, • • • ,d is the ith component of the parameterization function $ and the parame¬ 


terization function $ = <bj( y ), J(y) is the index function given in (3.12). a = a(x, y) = <hj^(y), 

£ = £(x, y) = (x) — $7 ( y)(y)- I R rest of the proof, without introducing any confusion, we 

always to use these short notations to save the space. In above derivation, we need the convexity 
property of the parameterization function to make sure all the integrals are well defined. 

Using above equality and the smoothness of the parameterization functions, it is easy to show 
that 


'Oi yMl, 


« i (x,y)|d(x,y)| 2 d^ x J d/z y 


< Ct 3 

< Ct 3 max 



0 JO JO JOi J My 


Rt(x,y) \D 2 ’ 3 u(<5> J(y) (a + s 3 S 2 Si^))| 2 d/i x d/i y ds 3 ds 2 dsi 
Rt(x, y) I D 2 ' 3 u($i(a + s£))| 2 d/i x d/r y , 

where we use the fact that J(y) = i, y G Oi and 


0<s<l J 0i J My 


|ZJ 2 ’ 3 u(x)|“ = lV j "V j 'V j u(x)l 2 + |V J "V J «(x)| 2 . 

j ,j',j"=l j,j'=i 

Let z i = <h ? ;(a + s£), 0 < s < 1, then for any y G Oi C B^. and x G _A4 y , 

|z* - y| < 2s|£| < 4s|x - y| < 8sVt, |z, - q*| < |z* - y| + |y - q*| < 5 + 8 sVt. 
We can assume that t is small enough such that 8 y/t < 5, then we have 


z i G B. 


25 

Ui' 
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After changing of variable, we obtain 


>Oi JM\ 


Rt{x, y) | D 2 ' 3 u($i(a + sO)| 2 dMxd/[ij 

f 1 „{\ Z i~y\ 2 \ ln2,3„./„ .m2 


< <± [ [ 1 p/ >i-y|' 

<5o JOi Jb z s s k V 128s 2 t 


D 2 , 3 u(zi) | d/i Zi d/i y 


2 / 


do An. \ 128s 2 f 


d/Jy / |-D 2,3 n(zj) I" d/x Zi 


< C f |A> 2 ’ 3 m(x)|^ d^x. 


This estimate would give us that 


|n(x)|| L 2 ( _ M) < Cf 1/2 ||w||^3 (A ^ ) 


(3.14) 


Now, we turn to estimate the gradient of ri. 


2 

[ |V x ri(x)| 2 d/r x < C f f V x Rt(x,y)d(x,y)p(y)dp y d/x a 
ijK JM JM 


!M {JM 


Rt(x, y)V x d(x, y)p(y)d/j y d/i„ 


where V x is the gradient in Ai with respect to x. 

Using the same techniques in the calculation of ||ri(x)|| L 2 (_ A/( ), we get that the first term of right 
hand side can bounded as follows 


lM \JM 


V x Rt(x, y)d(x, y)p(y)dp y dp x < C\[u\\ 2 h3{M) . 


The estimation of second term is a little involved. First, we have 


/ [ Rt(x, y)V x d(x, y)p(y)dp y dp x < C [ ([ R t (x,y)\V x d(x,y)\ 2 dpy] dp x 

JM JM JM \JM J 

= C'V [ ([ R t (x, y)\V X d{x, y)\ 2 dp^\ dp y 

i=l J °i V M y J 


Also using Newton-Leibniz formula, we have 


d(x, y) = CC J J si [di&(a + (a + s 2 si£)V j ' V j u($(a + s 2 si£))J ds 2 dsi 

-££' j J o s 1 (d i &(a)d i '&'{a)V j 'V j u{<!>(a))^ds 2 ds 1 
Then the gradient of d(x, y) has following representation, 


v x d(x, y ) = rrv. 


•Si (di&(a + si£)di>& (a + s 2 si£)V j ' V%(<f>(a + s 2 si£))) ds 2 dsi 


+V X (rr') J J J si-j^(di$ j {a + S3S 1 t,)di’$ j '(a + S3S2Si{,)'V j 'V j u($(a + s 3 S2Sit))'jds 3 ds2dsi 


= di(x,y) + d 2 (x,y). 


14 



For d\, we have 


[if i^(x,y)|di(x,y)| 2 d// x ) d/r y < Ct 2 max [if Rt(x, y)\D 2 ’ 3 u($i(a + s£))| 2 d/r x ) d/x y , 
JOi J 0<s<lJ o . yj M t y J 


which means that 


[([ J R f (x,y)|di(x,y)| 2 d^ x ) d^y < C / \D 2 ' 3 u(x)\ 2 dfjt* 

JOi \JM* y J Jb% 


(3.15) 


For d‘ 2 , we have 
^2(x,y) 

= V 


J 3 Sl ~^~ (^^(a + s 3 si£)di'$ J, (a + + S 3 S 2 S 1 O)) ds 3 ds 2 dsi 

= V x (Vf') C*" [ s\s 2 di&(a + S 3 Si€)di»di'$ j '(a + s 3 S 2 Si^)V i 'V :; u($(a + s 3 s 2 si£))ds 3 ds 2 dsi 

v y J[0,1] 3 

+V X (CC') C" [ s\di"di&(a + s 3 si£)di>&'(a + s 3 s 2 s 1 £'V j ''V j u(®(a + s 3 s 2 si£))ds 3 ds 2 dsi 

v ' J[0,1] 3 

+V X (CC') C" [ s?s 2 9j$ y (a + s 2 si^)5i/$ y, (a + S3S 2 si^)9j//$- 7 "(a + s 3 S2Si^) 

V ' J[0,1] 3 


V y V J V y tt(<I > (a + s 3 s 2 si^))ds 3 ds 2 dsi 


This formula tells us that 


[if i?t(x,y)|d 2 (x,y)| 2 d/x x ) d/r y < Ct 2 max [if #t(x,y)|T> 2 ’ 3 u($(a + s£))| 2 d// x ] d/x* 

JOi / o< s <i,y 0 . yv*,*, y 

Using the same arguments as that in the calculation of ||n.||z,2(.M)j we have 

[if i? i (x,y)|d 2 (x,y)| 2 d / u x ) d/x y < C [ |H 3 u(x)| 2 d/i x 

JOi \Jm $, / ^ 


Combining (3.15) and (3.16), we have 

l|Vn(x)|| L 2 ( _ M) < C\\u\\ H 3 {m) 

For r 2 , first, notice that 

V^(x, y) = i^^'(a)^ m '"X^ i (a)(* < - 2 / i )i 2 t(x,y), 

^t(x,y) = ^ m ^ y (a) 3 mV ^^(«)f'^^^(x,y). 


(3.16) 


(3.17) 


Then, we have 


V J i? t (x,y) - d_i? t (x,y) 


= 7 ^dm'$ l g m n dn'& ( x 3 - y 3 - £ d v & ) Rt(x, y) 


= Y/^ dm '^ 9m ' n ' dn ^ j {Jo Jo ^ 


sdjidi'&{a + rs£)drds ) i?*(x, y) 
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Thus, we get 


V J i?t(x,y) - ^i?t(x,y) < —^--Rt(x,y) 

V x fv J E t (x ) y)-^-i2 t (x,y) > ) < ®i?t(x,y) + ^S-|i^(x,y)| 


Then, we have following bound for r 2 , 


|r 2 (x)| 2 d/r x 


(3.18) 


M \J M 


R t (x,y)\D u(y)\p(y)dny ) d/i 


<Cf f ([ R t (x,y)p{y)dfj,y) [ Rt(x-,y)\D 2 u(y)\ 2 p(y)dpydpx 
JM \JM J JM 

<Ct,max(^J ^(x,y)d, x ) / |H 2 u(y)| 2 p(y)d / [i y 


<C7i 1111 ^-2 (jVi) - 
Similarly, we have 

f |Vr 2 (x)| 2 d/i x 


(3.19) 


X \JM 


V x i?i(x,y)p(y)d/r y ) [ V^R t (x,y)\D 2 u(y)\ 2 p(y)dp y dp^ 
J JM 


<CVtmax(^j V x i?t(x,y)d/i x ^ J \D 2 u(y)\ 2 p(y)dp y 

r 3 is relatively easy to estimate by using the well known Gauss formula. 

r 3 (x) = [ n J ?f(V*V J u(y))^ t (x, y)p(y)dr y - [ rf('V l V j u(y))R t (x, y)V j p(y)dp y 
JdM JM 

= hd ~ [ r/ ? (V*V J u(y))^ i (x,y)V J p(y)d / u y 
JM 

where I bd = J QM n-V(V*VMy))i?*(x, y)p(y)dr y . 

Using the assumption that p G C' 1 (A1), it is easy to get that 

ll r 3 - hd\\L 2 (M) < C , V / t||'u|| H 2 (_ M ), 

||V(r 3 - /m)||l 2 (>i) < C'||u||^2(_ M ). 

Now, we turn to bound the last term r 4 . Notice that 

V^'(VMy)) = (^^)/ r 9 r ((9 m ^> mV (a„rn)) 

= ( 4 '^)/ Z ' (fy( 3 m ^'))s mV ( 0 n'u) 

+{dk'&)g k ' l \dm'&)dv (g m ' n \d n 'uj) 


TtetG 

1 

/det G 


(3.20) 

(3.21) 

(3.22) 


(dm'V det G)g m ' n ' ( d n >u ) + <9 m ' ^ m ' n '(d n /u)) 
dm' (Vdet Gg m ' n '(dn'uf) = A M u(y). 
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where detG is the determinant of G and G = {g%j)i,j=i, - ,k- Here we use the fact that 

(d k ,&)g k ' 1 ' {d v (d m ^)) = (d k '&)g k ' 1 ' {d ml (d v ^)) 

= {d m ,{d k ,&))g k ' l \d v &) 

= ^g kl d m ’(g k 'i') 


-^(^VdetG). 


Moreover, we have 


g* i \d j .&)(d i 'Z l )(d l &)(y i V*u(y)) 

= -/ / (^'^)(^<h i )(V*V%( y)) 

= -g*1\d jl &)(d i '&)(d m '&)g m ' n 'd nl (V^(y)) 
= (V'u(y)) 

= -V J ' (V J u(y)) . 

where the first equalities are due to that d t ’l} = —5\,. Then we have 

div (r 7 *(V i V J u(y))) + X M u(y) 

1 


(3.23) 


VdetG 

? 

VdetG 


di> VdetG g'i (a,-^)e‘(^^)(V*V J «(y)) - ^ [d j ,&)[d i .g){d l &){V'V 3 u( y)) 


d v ( v det G g l 3 (a,-/^)(3i$*)(V‘V J tx 


Here we use the equalities (3.22), (3.23), 77 * = £ <9j/4r and the definition of div, 

divX = 1 fy(Vdet Gg i ' i 'd j '$ k X k ). 

V det G 


(3.24) 


where X is a smooth tangent vector field on A4 and (X 1 ,..., X d ) l is its representation in embedding 
coordinates. 

Hence, 


r 4 (x) = 


e 1 


dA VdetGg^ (d j ^)(di^)(V^u(y)))R t ^,y)p(y)dg. 


Im VdetG 

Then it is easy to get that 

ll r ’ 4 (x)|| L 2 ( _ M) < Ct 1/2 \\u\\ H 3 (M) , 
l|Vr 4 (x)|| i 2 (-M) < C\\u\\ H 3 {m) . 

By combining (3.14),(3.1T),(3.18),(3.19), ( |3.20[ ) , ( f3.21[ ),( |3.25[ ),( |3.26|) , we know that 

\\ r - hd\\L 2 (M) < Ct^ 2 \\u\\ H 3^ M ), 

IIV(r* — hd)\\L 2 (M) < c \\ u \\h 3 (M)- 
Using the definition of Im and Im, we obtain 


(3.25) 

(3.26) 


(3.27) 

(3.28) 


Ibd Ibd — 


IdM 


n J ( y)(x —y-7/(x,y)) • (VV J u(y))i? 4 (x,y)p(y)dr y 
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Using the definition of ry(x, y), it is easy to check that 


x —y — r/(x,y)| = 0(|x - y| 2 ), |V x (x - y - r/(x,y))| = 0(|x - y|) 


which implies that 


\\hd~ hd\\L 2 (M) < 

||V(/ M -/ M )|| L 2 {A1 ) < C't 1/1 ||li|| H 3(^ / (). 
The theorem is proved by putting ( 3.27[ >, ( |3.28[ ), ( |3.29[ ), ( |3.30[ ) together. 


(3.29) 

(3.30) 

10 


Remark 3.4. Using above proof, we can also show that the L 2 error in the integral approximation 
(2T3| is Oft 1/4 ). 


3.4 Proof of Theorem 13.31 

To simplify the notation, we introduce a intermediate operator defined as follows, 


Lt,hufx) 


7 R t(x,xj)(ufx) - u (xj))p(x j )V j . 

XjSP 


(3.31) 


Let ut t h = /f(u) with u satisfying equation (3.2) and If is given in (3.3). One can verify that the 
following equation are satisfied, 


-L t , h ut,h{x) = ^2 Rtfx,Xj)f(xj)/p{xj)Vj. 

x ,eP 


(3.32) 


In the proof, we need a prior estimate of u which is given as following. 


[pi - 

Theorem 3.6. Suppose u = ,U|p|) with Yji=\ u iPiVi = 0 solves the problem (3.2) and 

f = (/(x 1 ), • • • , /(x|p|))* for f € C(M). Then there exists a constant C > 0 such that 


(10 \ 1/2 

<c\\f\\oo 


provided t and h<yP '2/i M ^ are sma ^ enough. 

This theorem is an easy corollary of following theorem. 

Theorem 3.7. If the manifolds A4 is C°° , there exist constants C > 0, Cq > 0 independent on t 

so that for any u = (ui, • • • ,ii|p|)* £ with ^[=1 u iPiYi = 0 an d for any sufficient small t and 
h{P,V,M ) 

Vi 

IP) IPI 

y, -Rt(Xj,Xj)(uj - Uj) 2 piPjViVj > C( 1 - Coh(P, V,M) ^ y-v u 2 p ,y.' 

i,j= 1 V *=1 


The proof of this theorem is given in the supplementary material which is a small modification 
of the proof of Theorem 9.1 in 


We are now ready to prove Theorem 3.3 
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Proof. To simplify the notation, we denote h = h(P, V, M) and n = |P| and denote 


ut,h(x.) 


1 


y RtfcxjhjPjVj y R t( x ’ x j)fj v j/pj 

yxj-eP x jSP 


(3.33) 


where u = (u\, • • • with UiPiVi = 0 solves the problem (3.2), fj = f(x.j) and w t} hi x ) = 

eP x j)PjVj- For convenience, we set 


a u( x ) = — ^ x) Rti^^UjPjVj, 


x,eP 


c t,/i(x) = - * R t (x.,x. j )f(x. j )V j /p j , 

t ’ h ' ' XjSP 

and thus = a t ,h + ( ‘t,h- 

First we upper bound ||Lt(ri t)?l ) - Lt,h(wt,h)||z, 2 ^). For ct,h, we have 


(3.34) 

(3.35) 


K-C'tCt./i - L t , h c t ,h) (x)| 

1 

M 


t 


< I ^t,h ( x ) | 


1 

+ ~t 


-Rt(x,y)(c tj/l (x) - c tj/l (y))p(y)d// y - ^ i? t (x, x J )(c^ h (x) - c t:h (-x. j ))p j V j 

VLj&P 

^t(x,y)p(y)d^ y - -Rt(x, x j )p j V j 

x jSP 

i?t(x,y)c tj / l (y)p(y)d^ y - ^ R t (x.,x. :j )ct,h('Xj)PjVj 




Ch Ch 

— f3/2 I c t,/i( x )I T ^ 3/2 ll c *,^llc' 1 (A / () 

< %mu + + t 1/2 imioo) < noo. 


For at,h, we have 


< 


< 


< 


JM 

Ch 2 
t 


(at,h(x.)Y 




i?i(x,y)p(y)d^ y - y Rt(x,*j)PjVj 


d/r* 


(3.36) 


[ ( a t,h( x )) 2 d/r x < 

JM 


Ch 2 _ 

1 Jm \ ™t,h(x) 


^ Rtf*, Xj)ujPjVj dp* 



y />*/ (x. x y) iij pj V j y R. t (x, Xj )pj Vj d t iy 


Ch 2 


i2t(x,Xj)d/r x < —j—yUjPjVj. 


j= 1 
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Let 


A = C t [ R 

Jm w t,h{y) 


|X - y| 
4 1 


R 


IXi - y| 
At 


x,eP 


1 


i? 


wt,h(*j) V 4t 


x - X,- 


i? 


4* 


p(y) d h> 


Pj^- 


We have |A| < for some constant C independent of t. In addition, notice that only when 
|x - x, : | 2 < 16t is A ^ 0, which implies 


14 < 


|X - Xj| 

32t 


Then we have 






Rt(x,y)a t ,h{y)p(y)diiy - E J R t (x,x i )a tj / l (x i )p i V ? - 

XjSP 

2 


d hx 


^ \*=i 


Y CtUiPiVjA dp : 


< 




c_y 

t 


\ i=l 


^ ^ Ct\ui\piViR 


< 


Ch 2 
t 


j Y* 

Jm Vi.i 


<—t (/..«* 


t 


i=1 


'X 


|X ~ X, 

32t 


|X - X, 

32t 


u: 


l ViVi ) ( E 

XjGP 


X — X, 

32/ 


|X - Xj 

32t 
2 


d/J> 


Pi Vi dp. 


Ch 2 


dpx {ufpiVj) 1 < —— E n iP^ 


V Z=1 


Combining Equation (3.36), ( |3.37[ ) and Lemma 3.6 

ll-hiOpfo R t,h^t,h || L 2 (JA) 




I (Lt{a t ,h) ~ L t ,h(a t , h )) (x)| 2 dp x 


1/2 


1 

< - 
t 




(a t ,h(*)Y 


' M 


Rt{*, y)p(y)dp y - E /'’/(x. x/i/pi} 


x,eP 


2 X 1/2 

dp. 


(3.37) 


1 

+ t 


/X 


/X 


i7 i (x,y)a tj/l (y)p(y)dp y - E R Mi*j) a t,h(*j)Pj V j 


1/2 


dp. 


< 


Ch 


E u iP* y * 


1/2 


< 


^ 2=1 


Ch 


Assembling the parts together, we have the following upper bound. 

Lt : h^t,h\\L 2 (A4) 

Lt,h(H,h || L 2 (M) T \\ R tCt,h R t,hCt,h || L 2 (A4) 

Ch 


(3.38) 


Ch „ Ch 

< ^ll/lloo + — 


< 


*3/2 
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At the same time, since ut respectively ut.h solves equation (3.5) respectively equation (3.32), we 
have 


Lt{ut) Lf ! h[ut : h)\\L 2 (M) (3.39) 



< -^WfWcHM)- 


The complete L 2 estimate follows from Equation (3.38) and ( |3.39[ ). 

The estimate of the gradient, \\S7{L t {ut) — Tt,h{ut,h))\\L 2 {M)i can be obtained similarly. 


□ 


3.5 Proof of Theorem 13.21 


In order to prove Theorem |3.2[ we need two theorems, |3.8|and 3.9 The proof of these two theorems 


can be obtained by making minor revision of the proof of Theorem 4.4 and 4.5 in m, the details 
of the proof are put in the supplementary material. 


Theorem 3.8. For any function u £ L 2 (Ai), there exists a constant C > 0 independent on t and 
u, such that 


l i R ( ^ ) M( x ) ~ ^(y)) 2 p( x My)dh-xd/v > c [ |Vu(x)|Vx)d^ x , 

J M JM V 41 / JM 

where 

M x ) = C ! . [ R ( ^ ) u(y)p(y)dp. y , 

Wt(x) Jm V / 


and w t (x) = C t J M R p(y) d hy 

Theorem 3.9. Assume both AA. and dAi are C°°. There exists a constant C > 0 independent on 
t so that for any function u £ L,2(A4) with f M u(x)p(x)d/r x = 0 and for any sufficient small t 

J m J m r ) (M x ) _ ^(y)) 2 p( x )p(y) d Mxd^ y > C\\u\\l 2(M )- 


Using above two theorems, Theorem 3.2 becomes an easy corollary. 


Proof, of Theorem 3.2 

Using Theorem 3.9 we have 

12 


u \\l 2 {M) < C / u(x)r(x)p(x)dp x < C\\u\\ L 2 {M) \\r\\ L 2 {M) . 
J M 


This inequality (3.40) implies that 


Ml L 2 {M) < C'lkll L 2 (M)- 


(3.40) 
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Now we turn to estimate || Vti|| x, 2 (A^) - Notice that we have the following expression for u, 


1 i(x) = u(x) + 


^t(x) 


r x 


where 


v ( x ) = ~T^ f R t(x,y)u(y)p(y)dVy, w t (x) = R t (x,y)p{y)dn y . 

mW J m Jm 


By Theorem 3.8 we have 


ll Vn ll L 2 (M) - 2 II^ U II L 2 {M) + 


2 tr 


V 


r x — r 


wt(x) 


L 2 (M) 

2 | /o«2|i v7^I|2 


< C / u(x)L t u(x)p(x)d/r x + Ct||r|| jL 2 (A/f) ||Vr|| i2(>() 

J M 

< ^IMIl 2 (.M)IMI.L 2 (.M) + C^IMIl 2 (.M) + 

< C'||r||i a(A 4 ) +C't 2 ||Vr||i a{A<) 

< C (ll r llL 2 (X) + ^l|N7 r ||x 2 (7W)) • 


The proof is completed. 


3.6 Proof of Theorem 13.51 

Proof. First, we denote 


r(x) = / b(y) • (x - y )R t (x, y)p(y)dr y , 

JdM 

b(y) • (x - y)^t(x,y)p(y)dr y ) p(x)dx. 


r = 


IdM 

1 

\M\ 


P J M \JdM 


where \M\ P = f M p(y)dp y . 

The key point of the proof is to show that 




u(x) (r(x) — r)p(x)d/Xj 


< CVt \\b\\ H i^ M )\\u\\ H i( M y 


First, notice that 


|r| < CVt \\b\\ L 2 {dM) < CVt ||b||xri(7W)- 
Then it is sufficient to show that 


/ ^( x ) / b (y) • ( x - y)^t( x >y)F(y) dT y ) p( x ) d h> 

Im \JdM 


Direct calculation gives that 

|2tVi?t(x, y) - (x — y)Rf(x, y)| < C|x - y\ 2 R t (x, y), 


□ 


(3.41) 


< Cy/t ||b|| ii -i(^ 1 )||'u|| ii -i(^ 1 ). (3.42) 
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where Rt(x, y) = CtR j and R(r) = / r °° i?(s)ds. This implies that 

[ u(x)p(x) [ b(y) ((x - y)Rt(x,y) + 2tX7R t (x,y)) p(y)dr y d/x> 
JM JdM v 7 

<C [ \u{x)p(x)\ [ |b(y)||x - y| 2 J R t (x,y)p(y)dr y d/i x 


(3.43) 


>M 


<Ct\\h\\ L 2 {dM) ( I l^J ^R t {x,y)p(x)dp x J ^ j \u(x)\ R t (x,y)p(x)dp x ) p(y)dr y 

. 1 / 2 

<Ct\\b\\ H i(M) ( j M x )lV x ) \J KfryMdTy) d/i x 


<C't 3 / 4 ||b|| jff i(^ 1 )||?x|| L 2(_ M ). 

On the other hand, using the Gauss integral formula, we have 


1/2 


/ u{x)p(x) b(y) ■VRt(x,y)p(y)dT y dp^ 
>M JdM 


(3.44) 


/ax 4x 


w( x )p( x )T x (b(y)) • Vi? t (x,y)p(y)d/i x dr y 


IdM JdM 


n(x) • T x (b(y))tt(x)i?t(x, y)p(x)p(y)dr x dr y 


IdM JM 


div x [n(x)p(x)T x (b(y))]i?t(x, y)p(y)d/i x dr y . 


Here T x _is the projection operator to the tangent space on x. To get the first equality, we use the fact 
that VRt(x , y) belongs to the tangent space on x, such that b(y) ■ VRt(x, y) = T x (b(y)) ■ VRt(x, y) 
and n(x) • T x (b(y)) = n(x) • b(y) where n(x) is the out normal of dM at x G dM. 

For the first term, we have 


IdM JdM 


n(x) • T x (b(y))«(x)i? f (x,y)p(x)p(y)dr x drj 


(3.45) 


IdM JdM 

<C'l|b||L 2 (ax) 


n(x) • b(y)u(x)i? t (x,y)p(x)p(y)dr x dr y 

u( x )|fii( x ,y)p(x;)dT x ) p(y)dT y 


dM \JdM 


2 \ V2 


< C 'll b llH 1 (X) / / ^t( x ,y)p( x ) dT x / \u(x)\ i?t(x,y)p(x)dr x p(y)dr y 

\JdM \JdM J \JdM 

<Ct 1 / 2 \\b\\ H i(M)\\ u \\L 2 (dM) < Ct 1,/2 ||b||^i(x)||'u|| H i(_ M ). 


1/2 


We can also bound the second term on the right hand side of ( |3.44[ ) . By using the assumption that 
M G C °°, we have 

|div x [u(x)p(x)T x (b(y))]| 

<|Vu(x)||T x (b(y))||p(x)| + |«(x)||div x [T x (b(y))]||p(x)| + |Vp(x)||u(x)T x (b(y))| 
<C(|Vu(x)| + |«(x)|)|b(y)| 


where the constant C depends on the curvature of the manifold M. 
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Then, we have 


'dM J M 


div x [u(x)T x (b(y))].R i (x,y)p(x)p(y)d/x x dT y 


(3.46) 


<C b(y)p(y) / (|Vu(x)| + |n(x)|)iii(x,y)p(x)d^ x dr y 


IdM 




< C ||b|| L2(aM) / (|Vu(x)| +|«(x)| )p(x) I / i?i(x,y)p(y)dr y ) d/x, 

>M \JdM 


1/2 


< Ct l/ 1 \\b\\ H i^ M )\\u\\ H i^ M y 


Then, the inequality (3.42) is obtained from (3.43), (3.44), (3.45) and (3.46). Now, using Theorem 


3.9, we have 


\\ u \\h(M) - C u(x)L t u(x)p(x)dp x < CVt \\b\\ H i( M) \\u\\ H i {M) . 
J M 

Note r(x) = f dM {x — y) • b(y)^(x,y)p(y)dr y . Direct calculation gives us that 

ll^( x )II l 2 (A 1 ) < Ct 1/l \\b\\ HHM) , and 
l|Vr(x)|| L 2 (A ^) < C , i _1/4 ||b|| H i (A1) . 

The integral equation Ltu = r — r gives that 


u(x) = v(x) + ^ (r(x) — r) 


where 


V ( X ) = —TC Rt(*,y)u{y)p(y)dp y , w t (x)= R t {x, y)p{y)dp y . 
w t\. x ) JM JM 


By Theorem 3.8 we have 

l! Vu lli 2 (M) 

< 2||Vn||| 2 (_ A/ () + 


2t/ 


V 


r x — r 


w t {x) 


L 2 (M) 

2 i /^+2 ||y7„||2 


< C / u(x)L t u(x)p(x)dp x + Ct\\r\\ L2{M) + Ct \\Vr\\ L2[M) 

J M 

< Cy/t ||b||#i(_y\4)||w||_ffi(y4) + Ct\\r\\\ 2{M) + Ct 2 \\X 7 r\\ 2 L2 ( M ^ 

< ^1|b||//i(x) (yt\\u\\m( M ) + Ct i ^ 2 ^j . 


Using ( |3.47[ ) and ( |3.48[ ), we have 

IMI 2 m(M) — C||b||iJi(At) ^V^||“||j?i(x) + Ci 3/2 ) > 


(3.47) 


(3.48) 


which proves the theorem. 


□ 
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4 Numerical Experiments 


In this section, we show several numerical examples to demonstrate the performance of the point 
integral method for isotropic elliptic equations. This section is separated to two parts. In the first 
part, on some simple 2D surfaces, the convergence of the point integral method is verified. In the 
second part, we consider a nonlocal total variation minimization problem, in which some isotropic 
elliptic equations are solved on point cloud in high dimensional space. 

4.1 Examples on 2D Surfaces 

In this subsection, we consider the isotropic elliptic equation on 2D surfaces 

—div(p 2 (x)Vri(x)) = /(x), x e M, (4.1) 

with Neumann and Dirichlet boundary conditions, 

— (x) = 6(x), or uCx) = 6(x), x E DM 

on 

To estimate the volume weight vector V from the point sets P, a local mesh around each sample 
point is constructed, from which the weight of that point is computed. For details to estimate the 
volume weight, we refer to [26] . The kernel function is chosen to be Gaussian function, 

ii ‘ (x ' y, = (4^ exp (‘ M if L )' 

The parameter t is set as t = X)i=i p( x i)) > where p(x*) is the radius of 10 nearest neighbors 
of x, . 

Example 1 In the first example, the manifold M is an unit disk and an annulus in M 2 . The inner 
radius of the annulus is 1 and outer radius is 3. The exact solution is set to be n 9 j(x) = cos(27r||x||) 
in unit disk and u g t = sin(x + y) in the annulus, see Figure [I] The coefficient of the equation in 



(4.1) is 


P = 



(4.2) 
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IPI 

684 

2610 

10191 

40296 

disk 

0.364597 

0.214960 

0.111961 

0.056028 

annulus 

0.036760 

0.012227 

0.005557 

0.003542 


Table 1: I 2 error for u g t = cos(27rr) in the unit disk and sin(x + y ) in the annulus. 


both in the unit disk and annulus. The Neumann boundary condition is enforced in unit disk and 
we consider the Dirichlet boundary condition in the annulus. 

Table [l] list the I 2 error of the point integral method as the number of points grows. This result 
clearly shows the convergence of the point integral method. The convergence rate in I 2 error is 
approximately 1 / \/\P\- 

The eigenvalue problem with homogeneous Neumann boundary condition is also solved in the 
annulus. 


—div(p 2 (x)Vu(x)) =Au(x), 


du 

<9n 


(x) =0, 


x e M 
x G dM 


and the coefficient p is given in (4.2). 

The first 20 eigenvalues are plotted in Figure [2] The eigenvalues given by finite element method 
in the finest mesh is used as the true solution. Our result shows that the eigenvalue computed in 
the point integral method also converge. 



Figure 2: First 20 eigenvalues in the annulus with Neumann boundary condition with different 
point cloud. 


Example 2 Now, we solve equation (4.1) with Neumann condition and Dirichlet condition on a 
curved surface in M 3 . Let A4 be a cap on the unit sphere, whose height is 1/2 and the cap angle is 
7t/3, as shown in Figure [3j The coefficient of the equation is also given in (4.2). 

We set the ground truth to be u g t = x + y + z, where ( x,y, z ) is the coordinate in M 3 . 

The I 2 errors of the point integral method are listed in Table [2] The convergence rate for both 
boundary value problems are 1. 
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1.5 



Figure 3: Ground truth: u g t = x + y + z on the cap. 


\p\ 

1199 

4689 

18540 

73757 

Neumann 

0.036779 

0.015355 

0.007479 

0.003189 

Dirichlet 

0.007238 

0.001921 

0.001278 

0.000750 


Table 2: I 2 error for u g t = x + y + z on cap. 


The first 20 eigenvalues are also computed for homogeneous Neumann condition as shown in 
Figure [4} As the number of points increases, the eigenvalues given by PIM converge to those 
computed by FEM, which suggests the convergence of the point integral method. 



Figure 4: First 20 eigenvalues on the cap with Neumann boundary condition. 
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Example 3 In this example, we consider a more complex surface, a human face called ’’Alex”. 
The surface is sampled by 10597 points (Figure[5]) and the analytic form of the surface is not known. 
The coefficient of the equation in (|4.1[) is 


P 


2 


1 

sin(r/10)/2 + 1 


where r = \Jx 2, + y 2 + z 2 . In this example, we solve the eigenvalue problem of the isotropic elliptic 



Figure 5: (a) Face of Alex; (b) Coefficient: restriction of p 2 on A4 

operator. Several eigenfunctions computed by the point integral method are shown in Figure [6} 
From the examples in 2D surfaces, we see that PIM solves isotropic elliptic equations with 
Neumann and Dirichlet boundary very well. Moreover, the convergence rate is higher than that 
obtained in the convergence analysis. The point integral method is applicable to point cloud in high 
dimensional space, not only on the 2D surfaces. Next, we will show a high dimensional example. 

4.2 Nonlocal Total Variation Extension 

In this example, we consider an L\ extension on point cloud. The point cloud is constructed 
by using the patches of a 512 x 512 image, which is shown in Figure ]7](a). The original image is 
subsampled and only retain 10% of the pixels at random. The subsampled image is shown in Figure 
0b). One classical problem in image processing is to recover the image from the subsampled image. 
Here, rather than give an image reconstruction method, we only use this example to demonstrate 
the performance of the point integral method for isotropic elliptic equations. 

In this example, the point cloud consists of the patches of the original image. For each pixel Xi 
in the image /, we extract a patch around it of size 5x5 which is denoted as p Xi (f ), where / is 
the original image. Totally, we get 512 2 patches and each patch is 5 x 5. The collection of all the 
patches give a point cloud in M 25 . Denote this point cloud as P = {p Xi (f ) : i = 1, • • ■ , 512 2 }. The 
image is actually corresponding a function u on the point cloud P with u(p Xi (f )) = f(xi), f(xi) 
is the value of image / at pixel X{. Corresponding to the subsampled image, the value of function 
u is only known in the patches around the sampled pixels. The collection of all these patches is 
denoted as S. 
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Figure 6: Eigenfunctions on ‘Alex’ with homogenous Neumann boundary condition. 


Recently, manifold model attracts many attentions in image processing m- In manifold model, 
the point cloud P is assumed to be a sample of an underlying manifold, which is called patch 
manifold. The total variation is used as a regularization to reconstruct the image. The main idea 
is to minimize the total variation in the patch manifold, i.e., 

min || Vu||£i(_A/j), subject to: rt(x) =/(x), x S S. (4-3) 


The variation approach tells us that the optimal solution of (4.3) is given by solving following PDE, 

div ( j = o. 

viv«(x)i; 


with the Dirichlet type boundary condition 

«( x ) = /( x ), x e S. 

One natural method to solve above PDE is an iterative scheme, 

diV (WRr) = °’ ^ n+1 ( x ) = /( x ), xeS. (4.4) 

In each step, we need to solve an isotropic elliptic equation. 

Here, the gradient is computed by using an integral approximation also. 


V'u(x) 


1 

tw t (x.) 



Rt(*, y)(x - 


y)M x ) 


n(y))d^ y 
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u)t(x) = f M Rt{x, y)d/i y . In the computation, to avoid degenerate of the ellipticity, we regularize 
the coefficient by adding a small constant in the denominator, i.e., replace |Vu n (x)| by |Vu n (x)| +e 
in (4.4) with e = 10 -3 . The point cloud is assumed to be uniformly distributed, so the volume 
weight is uniform. The kernel function is Gaussian function. In this example, we use the integral 
approximation (2.13) with adaptive f(xj) = p(xj) 2 , where p(xj) is the radius of 20 nearest neighbors 

Of Xj. 



Figure 7: (a): original data; (b): 10% subsampled data. 

Figure [8](|a) shows the image reconstructed by L\ extension and Figure [8^b) gives the difference 
between the original image, Figure [T^a) and the reconstructed image Figure [8](|a). As we can see, 
L i extension gives very good reconstruction. This result shows that the point integral method solve 
the isotropic elliptic equation very well on point cloud. 


5 Conclusion 

In this paper, we generalize the point integral method to solve the isotropic elliptic equation. The 
point integral method is very easy to implement on point cloud, since it only needs the point cloud 
without any extra information. Moreover, it also has very good theoretical property. The coercivity 
of the original elliptic operator is partially preserved in the point integral method. Based on this 
property, the convergence is proved. 

One important implication is the spectral convergence of the point integral method on random 
samples. Suppose the points are obtained by sampling a manifold according to some probability 
distribution p(x). In the point integral method, the eigenvalue problem 

j -^) div b 2 ( x ) Vu ( x )) = M x ), x€M, 

\ g*(x) = 0, xedM, [ ' j 
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Figure 8: (a): reconstructed data; (b): residual. 


is discretized as 


7 Rt(Xi,Xj)(Ui-Uj) 
x,eP 


A ^2 Rt(xi,*-j)uj- 


(5.2) 


This discretization is closely related with the normalized graph laplacian. Based on the theoretical 
results in this paper, it can be proved that the spectra of (5.2) converges to the spectra of (5.1) as 
the number of sample points goes to infinity. 

The other interesting problem is how to generalize the point integral method to anisotropic 
elliptic equation. On this problem, we already get some results. They are going to be reported in 
the subsequent paper. 
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